
clear all
addpath('..\functionsRA')

addpath('..\input')
addpath('..\input')

load('..\input\dists_US.mat','years','ny','last_month')
ending = '_US_monthquart';


load('..\input\results_est_101_US_monthquart.mat','Qs','pinew','R2')
Qs_q = Qs;

lth0 = size(Qs,1);

f=3;
lth = size(Qs_q,1);


% obs: in the graph above sometimes it seems that the quarterly is doing better, which is not possible, but that's just because those are months in which the quarterly model wasn't used (months 1,3,4,6,etc). 
% The dotted lines are the means of the fit in which both models are used, and monthly is slightly better. To check this, plot plot((1:44),[Qs((2:f:lth0-3),1) Qs_q]) and check they are basicaly the same, ans also that
% Qs((2:f:lth0-3),1) - Qs_q is alaways negative

mdl = 1;
name = ['model_' num2str(mdl+100)];

% graph with the prob of average inflation from t+6~t+10 being a disaster
pi_lim = [(-1:5)'; Inf]; pi_bar = [-2; (-0.5:4.5)'; 6];
prob_avg_dis = zeros(lth0,2); prob_avg_dis04 = prob_avg_dis; prob_avg_dis04_1to5 = prob_avg_dis; prob_avg_dis_1to5 = prob_avg_dis;
load('..\input\dists_US.mat','infl')
infl = infl(:);

%%
counter=1;
J=2;
simul_size=30*10^6;
Simuls=zeros(11,simul_size);
numSteps=10;
Cond_prob=zeros(4,32);
for i = 121:166
    Simuls=zeros(11,simul_size);

    bin = (pi_lim>=infl(i)); bt = find(bin,1);
    x = pinew{1,1}(:,i);
    A = Amatrix(x,mdl);
    mc = dtmc(A);
    x0=zeros(1,8);
    x0(1,bt)=simul_size;
    Simuls = simulate(mc,numSteps,'X0',x0);
    Check=(Simuls==bt);
  %% unconditional probability
    
All_steps=Simuls(7:11,:);
All_values=zeros(5,simul_size);
% 
old_values = 1:8;
new_values = pi_bar;
% 
for k = 1:length(old_values)
     All_values(All_steps == old_values(k)) = new_values(k);
end
Cond_prob(J+2,counter)=sum((log(prod(1+All_values/100))/5)>0.04)/simul_size;
% 
Cond_prob(J+1,counter)=sum(mean(All_values)>4)/simul_size;
Cond_prob(J+2,counter)=sum(mean(All_values)<0)/simul_size;

    %%
    for j=1:J
        % Which columns satisfy the condition
        good_columns=(sum(Check(2:j+1,:),1)==j);
        
        % probabbality of condition
        Prob_B=sum(good_columns)/simul_size;
      
        
        % Joint probability
        Future=Simuls(7:11,good_columns);
        Future_values=zeros(5,length(good_columns));
old_values = 1:8;
new_values = pi_bar;

for k = 1:length(old_values)
    Future_values(Future == old_values(k)) = new_values(k);
end
Check2= mean(Future_values)>4;
good_columns2=sum(Check2);
Prob_AB=sum(good_columns2)/simul_size;

Cond_prob(j,counter)=Prob_AB/Prob_B;
 end
    
    counter=counter+1;
 
end

currentFolder = pwd; 
subfolder = fullfile(fileparts(currentFolder), 'input'); 
mydatafinalUS = fullfile(subfolder,'myData_USA_FINAL.xlsx');

writematrix(Cond_prob', mydatafinalUS);
